clear all;
%close all;
clc;
%using a Matlab interpolation to get intermediate values for population
%using Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)
%and getting the predicted income series from 2001 onward 
%formally known as interpolation17

SN=50; %state numbers (here we do work without DC)
YEARS=42; %making number of years a variable
load ('input1960-2001all.mat');

for i=1:YEARS
    temp(:,i)=pop(SN*(i-1)+1:SN*(i));
end
clear pop
pop=temp;
clear temp;

for i=1:YEARS
    temp(:,i)=income(SN*(i-1)+1:SN*(i));
end
clear income
income=temp;
clear temp;

for i=1:YEARS
    temp(:,i)=incomepc(SN*(i-1)+1:SN*(i));
end
clear incomepc
incomepc=temp;
clear temp;

[pop04_30] = xlsread('pop2004_2030.xls');  
pop04_30=pop04_30./1000000;    % population in millions 


income(51,:)=sum(income);

for i=1:SN
x = [2000,2001,[2004:2030]]; 
y(i,:) = [pop(i,41),pop(i,42),pop04_30(i,:)]; 
t = [2000:2030];
p(i,:) = pchip(x,y(i,:),t);
end;


popall=[pop(:,1:40) p];
popall(51,:)=sum(popall(1:SN,:));


[area]=xlsread('area.xls');  %are in square miles
Area=kron(ones(1,71),area);
popdens=(popall(1:51,:)*1000000)./Area;
Popdens=popdens;

Popall{1}=popall;
Popall{2}=popall;
Popall{3}=popall;


%growth rates for population in sample and out of sample 
            
   for i=2:71
      popgth(:,i-1)=(popall(:,i)-popall(:,i-1))./popall(:,i);
   
  end 
       clear i;
       

%growth rates of income in sample

 for i=2:42
      incomegth(:,i-1)=(income(:,i)-income(:,(i-1)))./income(:,i);
  end 
       clear i;
 popgrowth={popgth popgth popgth};   
            
% Create bivariate normal distribution for Population and GDP Growth Rate          
for q=1:3
popgrth=popgrowth{q};
popall=Popall{q};
for i=1:SN   
    pm=mean(popgth(i,1:41));
    ps = std(popgth(i,1:41));
    
 popstd={ps ps ps};   
 popmean={pm pm pm};  
    
    g=mean(incomegth(i,1:41));
    gs=std(incomegth(i,1:41));
incmean={(g-0.01) (g) (g+0.01)};
incstd={(gs-0.015)  gs (gs+0.015)};
    
    rho= corrcoef([popgth(i,1:41)', incomegth(i,1:41)']);
    
    rho=rho(1,2);

ps=popstd{q};
p=popmean{q};
g=incmean{q};
gs=incstd{q};
    
    
    beta=rho*ps*gs/(ps*ps);
            
            for v=1:29;
                incomegth(i,(41+v))=(g-beta*pm+beta*popgrowth{q}(i,(41+v))) + (gs*sqrt(1-(rho*rho)))*randn(1);
            end;
            
end
clear i v
   % Calculating income Projection for all states
            for i=1:29;
                income(1:SN,(42+i))=income(1:SN,(41+i)).*(ones(SN,1)+incomegth(1:SN,(41+i)));
            end;
            clear i;
            
  %aggregate income for ploting
  income(51,:)=sum(income(1:SN,:));
  for i=2:42
  incomegth(51,i-1)=(income(51,i)-income(51,(i-1)))./income(51,i);
    end
      g_us=mean(incomegth(51,1:41));
    gs_us=std(incomegth(51,1:41));  
    
  
  %calculating per capita income in each state          
    for i=1:29           
  incomepc(1:SN,(42+i))=income(1:SN,(42+i))./(popall(1:SN,(42+i)));
    end
clear i

incomepc(51,:)=income(51,:)./(popall(51,:));
Income{q}=income;
Incomepc{q}=incomepc;
end

